Fit density models to cod and flounder and check spatial overlap

Intro

In this script, I take the collated stomach data set and calculate aggregates (feeding ratio, total weight of prey groups) and predictor variables for diet data, aggregate to get 1 stomach = 1 row per prey type (not prey individual). I also select only the columns I need for model fitting, join environmental covariates and cpue covariates for cod and flounder, and lastly saduria biomass densities.

Load packages & source functions

# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(patchwork)
library(janitor)
library(forcats)
library(gapminder)
library(viridis)
library(ggridges)
library(raster)
library(icesDatras)
library(ggalluvial)
library(ggrepel)
library(ncdf4)
library(chron)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(geosphere)
library(quantreg)
#> Warning in .recacheSubclasses(def@className, def, env): undefined subclass
#> "numericVector" of class "Mnumeric"; definition not updated
library(brms)
#> Warning: package 'Rcpp' was built under R version 4.0.5
library(sdmTMB)
library(sp)
library(raster)
library(ggcorrplot)
library(ggh4x)
options(mc.cores = parallel::detectCores()) 

# Source code for map plots
source("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/R/functions/map_plot.R")

theme_set(theme_plot())

# Continuous colors
options(ggplot2.continuous.colour = "viridis")

# Discrete colors
scale_colour_discrete <- function(...) {
  scale_colour_brewer(palette = "Dark2")
}

scale_fill_discrete <- function(...) {
  scale_fill_brewer(palette = "Dark2")
}

Read and summarise data

cpue_full <- readr::read_csv("data/clean/catch_by_length_q1_q4.csv") %>%
  janitor::clean_names() %>% 
  rename(X = x,
         Y = y) 
#> Rows: 1759609 Columns: 22
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (8): Country, haul.id, IDx, ices_rect, id_haul_stomach, species, haul.i...
#> dbl (14): density, year, lat, lon, quarter, Month, sub_div, length_cm, depth...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> rename: renamed 2 variables (X, Y)

# Summaries density by species group an haul
cpue <- cpue_full %>%
  mutate(species2 = species,
         species2 = ifelse(species == "cod" & length_cm >= 25, "large_cod", species),
         species2 = ifelse(species == "cod" & length_cm < 25, "small_cod", species2)) %>% 
  group_by(haul_id, species2) %>% 
  summarise(density = sum(density)) %>% 
  ungroup()
#> mutate: new variable 'species2' (character) with 3 unique values and 0% NA
#> group_by: 2 grouping variables (haul_id, species2)
#> summarise: now 27,897 rows and 3 columns, one group variable remaining (haul_id)
#> ungroup: no grouping variables

# Trim data with quantiles
ggplot(cpue, aes(density)) + 
  geom_histogram() +
  facet_wrap(~species2, scales = "free")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

cpue <- cpue %>% 
  group_by(species2) %>% 
  mutate(dens_quants = quantile(density, probs = 0.999)) %>% 
  filter(density < dens_quants) %>% 
  ungroup() %>% 
  dplyr::select(-dens_quants)
#> group_by: one grouping variable (species2)
#> mutate (grouped): new variable 'dens_quants' (double) with 3 unique values and 0% NA
#> filter (grouped): removed 30 rows (<1%), 27,867 rows remaining
#> ungroup: no grouping variables

# Make wide
wcpue <- cpue %>% pivot_wider(names_from = species2, values_from = density)
#> pivot_wider: reorganized (species2, density) into (flounder, large_cod, small_cod) [was 27867x3, now 9373x4]

head(wcpue)
#> # A tibble: 6 × 4
#>   haul_id                  flounder large_cod small_cod
#>   <chr>                       <dbl>     <dbl>     <dbl>
#> 1 1993:1:GFR:SOL:H20:21:1     9.41      77.9    0.628  
#> 2 1993:1:GFR:SOL:H20:22:32    6.63       5.13   0      
#> 3 1993:1:GFR:SOL:H20:23:31    0.953      2.61   0.00459
#> 4 1993:1:GFR:SOL:H20:24:30    1.89       9.71   0      
#> 5 1993:1:GFR:SOL:H20:25:2    16.7      400.     4.78   
#> 6 1993:1:GFR:SOL:H20:26:3    13.5      179.     2.95

# Left join in trawl information
nrow(wcpue)
#> [1] 9373
hauls <- cpue_full %>% distinct(haul_id, .keep_all = TRUE) %>% dplyr::select(-density)
#> distinct: removed 1,750,236 rows (99%), 9,373 rows remaining
nrow(hauls)
#> [1] 9373

cpue2 <- left_join(wcpue, hauls)
#> Joining, by = "haul_id"
#> left_join: added 20 columns (year, lat, lon, quarter, country, …)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 9,373
#> > =======
#> > rows total 9,373

colnames(cpue2)
#>  [1] "haul_id"         "flounder"        "large_cod"       "small_cod"      
#>  [5] "year"            "lat"             "lon"             "quarter"        
#>  [9] "country"         "month"           "i_dx"            "ices_rect"      
#> [13] "sub_div"         "length_cm"       "id_haul_stomach" "species"        
#> [17] "haul_id_size"    "substrate"       "depth"           "temp"           
#> [21] "oxy"             "sal"             "X"               "Y"

cpue2 <- cpue2 %>% drop_na(oxy, temp, depth, sal, substrate, flounder, small_cod, large_cod)
#> drop_na: removed 405 rows (4%), 8,968 rows remaining

# Scale variables
cpue2 <- cpue2 %>% 
  mutate(depth_ct = depth - mean(depth),
         depth_sq = depth_ct * depth_ct,
         depth_sq_sc = (depth_sq - mean(depth_sq)) / sd(depth_sq),
         depth_sc = (depth - mean(depth)) / sd(depth),
         temp_ct = temp - mean(temp),
         temp_sq = temp_ct * temp_ct,
         temp_sq_sc = (temp_sq - mean(temp_sq)) / sd(temp_sq),
         temp_sc = (temp - mean(temp)) / sd(temp),
         oxy_sc = (oxy - mean(oxy)) / sd(oxy),
         sal_sc = (sal - mean(sal)) / sd(sal),
         fyear = as.factor(year),
         fsubstrate = as.factor(substrate))
#> mutate: new variable 'depth_ct' (double) with 149 unique values and 0% NA
#>         new variable 'depth_sq' (double) with 149 unique values and 0% NA
#>         new variable 'depth_sq_sc' (double) with 149 unique values and 0% NA
#>         new variable 'depth_sc' (double) with 149 unique values and 0% NA
#>         new variable 'temp_ct' (double) with 8,679 unique values and 0% NA
#>         new variable 'temp_sq' (double) with 8,679 unique values and 0% NA
#>         new variable 'temp_sq_sc' (double) with 8,679 unique values and 0% NA
#>         new variable 'temp_sc' (double) with 8,679 unique values and 0% NA
#>         new variable 'oxy_sc' (double) with 8,678 unique values and 0% NA
#>         new variable 'sal_sc' (double) with 8,690 unique values and 0% NA
#>         new variable 'fyear' (factor) with 28 unique values and 0% NA
#>         new variable 'fsubstrate' (factor) with 5 unique values and 0% NA

Read and scale the prediction grid

pred_grid_1 <- read_csv("data/clean/pred_grid_(1_2).csv")
#> Rows: 296520 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (2): substrate, ices_rect
#> dbl (16): X, Y, depth, year, lon, lat, quarter, oxy, temp, sal, sub_div, den...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid_2 <- read_csv("data/clean/pred_grid_(2_2).csv")
#> Rows: 296520 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (2): substrate, ices_rect
#> dbl (16): X, Y, depth, year, lon, lat, quarter, oxy, temp, sal, sub_div, den...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

pred_grid <- bind_rows(pred_grid_1, pred_grid_2)

# Scale with respect to data!
pred_grid <- pred_grid %>% 
  drop_na(substrate) %>% 
  mutate(X = X / 1000,
         Y = Y / 1000,
         depth_ct = depth - mean(cpue2$depth),
         depth_sq = depth_ct * depth_ct,
         depth_sq_sc = (depth_sq - mean(cpue2$depth_sq)) / sd(cpue2$depth_sq),
         depth_sc = (depth - mean(cpue2$depth)) / sd(cpue2$depth),
         temp_ct = temp - mean(cpue2$temp),
         temp_sq = temp_ct * temp_ct,
         temp_sq_sc = (temp_sq - mean(cpue2$temp_sq)) / sd(cpue2$temp_sq),
         temp_sc = (temp - mean(cpue2$temp)) / sd(cpue2$temp),
         oxy_sc = (oxy - mean(cpue2$oxy)) / sd(cpue2$oxy),
         sal_sc = (sal - mean(cpue2$sal)) / sd(cpue2$sal),
         fyear = as.factor(year),
         fsubstrate = as.factor(substrate))
#> drop_na: removed 280 rows (<1%), 592,760 rows remaining
#> mutate: changed 592,760 values (100%) of 'X' (0 new NA)
#>         changed 592,760 values (100%) of 'Y' (0 new NA)
#>         new variable 'depth_ct' (double) with 7,528 unique values and 0% NA
#>         new variable 'depth_sq' (double) with 7,528 unique values and 0% NA
#>         new variable 'depth_sq_sc' (double) with 7,528 unique values and 0% NA
#>         new variable 'depth_sc' (double) with 7,528 unique values and 0% NA
#>         new variable 'temp_ct' (double) with 592,704 unique values and 0% NA
#>         new variable 'temp_sq' (double) with 592,704 unique values and 0% NA
#>         new variable 'temp_sq_sc' (double) with 592,704 unique values and 0% NA
#>         new variable 'temp_sc' (double) with 592,704 unique values and 0% NA
#>         new variable 'oxy_sc' (double) with 592,704 unique values and 0% NA
#>         new variable 'sal_sc' (double) with 592,760 unique values and 0% NA
#>         new variable 'fyear' (factor) with 28 unique values and 0% NA
#>         new variable 'fsubstrate' (factor) with 5 unique values and 0% NA

# Split by quarter
pred_grid_q1 <- pred_grid %>% filter(quarter == 1)
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
pred_grid_q4 <- pred_grid %>% filter(quarter == 4)
#> filter: removed 296,380 rows (50%), 296,380 rows remaining

# Load models
# qwraps2::lazyload_cache_dir(path = "R/main_analysis/spatial_overlap_cache/html")

Plot data in space

# hist(cpue2$depth)
# hist(log(cpue2$depth))
# 
# cpue_long <- cpue2 %>%
#   pivot_longer(c(flounder, small_cod, large_cod)) %>% 
#   rename(species2 = name, density = value) %>% 
#   pivot_longer(c(depth, temp, oxy, sal)) %>% 
#   rename(env_var = name, env_var_value = value)
# 
# ggplot(cpue_long, aes(density)) + 
#   geom_histogram() + 
#   facet_wrap(~species2, scales = "free", ncol = 1)
# 
# ggplot(cpue_long, aes(env_var_value, density)) + 
#   geom_point(alpha = 0.3) + 
#   stat_smooth(method = "lm", formula = y ~ x, color = "blue", se = FALSE) + 
#   stat_smooth(method = "lm", formula = y ~ x + I(x^2), color = "red", se = FALSE) + 
#   facet_wrap(env_var~species2, scales = "free", ncol = 3) + 
#   coord_cartesian(expand = 0)
# 
# cpue2 <- cpue2 %>% mutate(fle_presence = ifelse(flounder == 0, "N", "Y"),
#                           l_cod_presence = ifelse(large_cod == 0, "N", "Y"),
#                           s_cod_presence = ifelse(small_cod == 0, "N", "Y"))
# 
# # Biomass density
# plot_map_fc + 
#   geom_point(data = cpue2, aes(X*1000, Y*1000, color = flounder)) + 
#   scale_color_viridis_c(trans = "sqrt") + 
#   facet_wrap(~year)
# 
# plot_map_fc + 
#   geom_point(data = cpue2, aes(X*1000, Y*1000, color = large_cod)) + 
#   scale_color_viridis_c(trans = "sqrt") + 
#   facet_wrap(~year)
# 
# plot_map_fc + 
#   geom_point(data = cpue2, aes(X*1000, Y*1000, color = small_cod)) + 
#   scale_color_viridis_c(trans = "sqrt") + 
#   facet_wrap(~year)
# 
# # Presence / absence
# plot_map_fc + 
#   geom_point(data = cpue2, aes(X*1000, Y*1000, color = fle_presence)) + 
#   facet_wrap(~year)
# 
# plot_map_fc + 
#   geom_point(data = cpue2, aes(X*1000, Y*1000, color = l_cod_presence)) + 
#   facet_wrap(~year)
# 
# plot_map_fc + 
#   geom_point(data = cpue2, aes(X*1000, Y*1000, color = s_cod_presence)) + 
#   facet_wrap(~year)

Create meshes

# Q1
spde_q1 <- make_mesh(filter(cpue2, quarter == 1),
                     xy_cols = c("X", "Y"),
                     n_knots = 150, 
                     type = "kmeans",
                     seed = 42)
#> filter: removed 3,638 rows (41%), 5,330 rows remaining

# All flounder Q4
spde_q4 <- make_mesh(filter(cpue2, quarter == 4),
                     xy_cols = c("X", "Y"),
                     n_knots = 150, 
                     type = "kmeans",
                     seed = 42)
#> filter: removed 5,330 rows (59%), 3,638 rows remaining

Fit models

q1

# Small cod model q1 spatial
mcod_s_q1_space <- sdmTMB(small_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                            temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                          data = filter(cpue2, quarter == 1),
                          mesh = spde_q1,
                          family = tweedie(link = "log"),
                          spatiotemporal = "IID",
                          spatial = "on",
                          time = "year",
                          reml = FALSE,
                          control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 3,595 rows (41%), 5,271 rows remaining

sanity(mcod_s_q1_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_s_q1_space, conf.int = TRUE)
#>                             term     estimate  std.error    conf.low
#> 1              fsubstratebedrock -0.691928335 1.02212747 -2.69526137
#> 2            fsubstratehard clay -0.132212753 0.63380026 -1.37443843
#> 3  fsubstratehard-bottom complex -0.125952616 0.65479459 -1.40932644
#> 4                  fsubstratemud -0.140027403 0.63378872 -1.38223047
#> 5                 fsubstratesand -0.294006302 0.63514592 -1.53886944
#> 6                      fyear1994 -0.189885126 0.61380222 -1.39291537
#> 7                      fyear1995  0.298079563 0.60713833 -0.89188969
#> 8                      fyear1996  0.314750404 0.61276989 -0.88625652
#> 9                      fyear1997 -0.669099743 0.61556223 -1.87557955
#> 10                     fyear1998 -0.040065644 0.60397679 -1.22383840
#> 11                     fyear1999 -0.248289902 0.59798496 -1.42031889
#> 12                     fyear2000  0.687378415 0.62234748 -0.53240024
#> 13                     fyear2001  0.455272793 0.59180912 -0.70465176
#> 14                     fyear2002  1.962095060 0.59529336  0.79534152
#> 15                     fyear2003  0.320007017 0.61294634 -0.88134573
#> 16                     fyear2004 -0.220180305 0.58766498 -1.37198251
#> 17                     fyear2005  1.788929063 0.56681351  0.67799500
#> 18                     fyear2006  0.259180987 0.58035478 -0.87829347
#> 19                     fyear2007  0.957228921 0.56744971 -0.15495207
#> 20                     fyear2008  1.371939968 0.56516783  0.26423137
#> 21                     fyear2009  0.976393169 0.56340389 -0.12785817
#> 22                     fyear2010  1.472598180 0.57147570  0.35252640
#> 23                     fyear2011  0.353731080 0.58291380 -0.78875897
#> 24                     fyear2012  0.323141901 0.57780572 -0.80933650
#> 25                     fyear2013  1.624300129 0.56495747  0.51700384
#> 26                     fyear2014  1.075135510 0.57314313 -0.04820438
#> 27                     fyear2015 -0.082199342 0.57483185 -1.20884907
#> 28                     fyear2016 -0.073825753 0.57662362 -1.20398728
#> 29                     fyear2017 -0.127635067 0.58130161 -1.26696528
#> 30                     fyear2018  0.629950035 0.57463000 -0.49630407
#> 31                     fyear2019 -0.031397966 0.62355142 -1.25353629
#> 32                     fyear2020 -1.173909468 0.62458372 -2.39807106
#> 33                      depth_sc  1.162865195 0.08949681  0.98745466
#> 34                   depth_sq_sc -1.222660927 0.07417559 -1.36804242
#> 35                       temp_sc  0.943382195 0.13622800  0.67638023
#> 36                    temp_sq_sc -0.472594618 0.07282391 -0.61532686
#> 37                        oxy_sc  0.561042569 0.08513249  0.39418595
#> 38                        sal_sc -0.001169043 0.10988511 -0.21653990
#>      conf.high
#> 1   1.31140470
#> 2   1.11001293
#> 3   1.15742121
#> 4   1.10217566
#> 5   0.95085683
#> 6   1.01314511
#> 7   1.48804882
#> 8   1.51575733
#> 9   0.53738006
#> 10  1.14370711
#> 11  0.92373909
#> 12  1.90715707
#> 13  1.61519734
#> 14  3.12884860
#> 15  1.52135976
#> 16  0.93162190
#> 17  2.89986313
#> 18  1.39665545
#> 19  2.06940992
#> 20  2.47964856
#> 21  2.08064451
#> 22  2.59266996
#> 23  1.49622112
#> 24  1.45562030
#> 25  2.73159642
#> 26  2.19847540
#> 27  1.04445039
#> 28  1.05633577
#> 29  1.01169515
#> 30  1.75620414
#> 31  1.19074036
#> 32  0.05025212
#> 33  1.33827573
#> 34 -1.07727944
#> 35  1.21038416
#> 36 -0.32986238
#> 37  0.72789919
#> 38  0.21420182

# Large cod model q1 spatial
mcod_l_q1_space <- sdmTMB(large_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                            temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                          data = filter(cpue2, quarter == 1),
                          mesh = spde_q1,
                          family = tweedie(link = "log"),
                          spatiotemporal = "IID",
                          spatial = "on",
                          time = "year",
                          reml = FALSE,
                          control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 3,595 rows (41%), 5,271 rows remaining

sanity(mcod_l_q1_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_l_q1_space, conf.int = TRUE)
#>                             term    estimate  std.error    conf.low   conf.high
#> 1              fsubstratebedrock  3.35027430 0.73692195  1.90593382  4.79461479
#> 2            fsubstratehard clay  3.86345204 0.47295798  2.93647144  4.79043264
#> 3  fsubstratehard-bottom complex  3.68218901 0.48661209  2.72844685  4.63593118
#> 4                  fsubstratemud  3.87755397 0.47296062  2.95056819  4.80453976
#> 5                 fsubstratesand  3.89542852 0.47391863  2.96656508  4.82429196
#> 6                      fyear1994  0.42180846 0.51298108 -0.58361598  1.42723290
#> 7                      fyear1995  0.17395551 0.51194085 -0.82943012  1.17734113
#> 8                      fyear1996  1.05226053 0.51303144  0.04673738  2.05778367
#> 9                      fyear1997 -0.43001992 0.51539334 -1.44017230  0.58013247
#> 10                     fyear1998 -0.58814899 0.51441941 -1.59639250  0.42009452
#> 11                     fyear1999 -0.73612259 0.51068879 -1.73705423  0.26480904
#> 12                     fyear2000 -0.44805678 0.53791590 -1.50235257  0.60623901
#> 13                     fyear2001 -0.05835698 0.51165985 -1.06119185  0.94447790
#> 14                     fyear2002  0.32081027 0.51994811 -0.69826929  1.33988984
#> 15                     fyear2003  0.29559541 0.52581232 -0.73497780  1.32616861
#> 16                     fyear2004 -1.54152654 0.51295582 -2.54690147 -0.53615161
#> 17                     fyear2005 -0.23565112 0.49440049 -1.20465827  0.73335602
#> 18                     fyear2006  0.42869636 0.49515684 -0.54179321  1.39918593
#> 19                     fyear2007 -0.22043688 0.49335112 -1.18738732  0.74651355
#> 20                     fyear2008  0.17717789 0.49054227 -0.78426729  1.13862306
#> 21                     fyear2009  0.44018278 0.48883021 -0.51790683  1.39827239
#> 22                     fyear2010  0.35962700 0.49485452 -0.61027005  1.32952404
#> 23                     fyear2011  0.16164975 0.49879137 -0.81596338  1.13926288
#> 24                     fyear2012 -0.45179824 0.50054518 -1.43284877  0.52925229
#> 25                     fyear2013 -0.06603709 0.49115060 -1.02867458  0.89660039
#> 26                     fyear2014 -0.18103119 0.49683537 -1.15481063  0.79274826
#> 27                     fyear2015 -0.13023130 0.49483422 -1.10008855  0.83962595
#> 28                     fyear2016 -0.20703546 0.49514267 -1.17749726  0.76342635
#> 29                     fyear2017 -0.55937561 0.49580712 -1.53113971  0.41238849
#> 30                     fyear2018 -0.46754421 0.49737326 -1.44237789  0.50728948
#> 31                     fyear2019 -1.04469284 0.53930616 -2.10171350  0.01232781
#> 32                     fyear2020 -1.30245419 0.52725884 -2.33586253 -0.26904585
#> 33                      depth_sc  0.75692314 0.06450849  0.63048883  0.88335745
#> 34                   depth_sq_sc -0.34635027 0.03128711 -0.40767189 -0.28502866
#> 35                       temp_sc  0.69470536 0.09776062  0.50309806  0.88631266
#> 36                    temp_sq_sc -0.18487760 0.05099597 -0.28482787 -0.08492734
#> 37                        oxy_sc  0.32515528 0.06476164  0.19822480  0.45208577
#> 38                        sal_sc -0.10165223 0.07887191 -0.25623833  0.05293388

# Flounder model q1 spatial
mfle_q1_space <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                          temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                        data = filter(cpue2, quarter == 1),
                        mesh = spde_q1,
                        family = tweedie(link = "log"),
                        spatiotemporal = "IID",
                        spatial = "on",
                        time = "year",
                        reml = FALSE,
                        control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 3,595 rows (41%), 5,271 rows remaining

sanity(mfle_q1_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfle_q1_space, conf.int = TRUE)
#>                             term     estimate  std.error      conf.low
#> 1              fsubstratebedrock  4.217305679 0.50885474  3.2199687197
#> 2            fsubstratehard clay  3.793488470 0.41832849  2.9735796929
#> 3  fsubstratehard-bottom complex  3.709849071 0.42664643  2.8736374301
#> 4                  fsubstratemud  3.924364749 0.41826303  3.1045842668
#> 5                 fsubstratesand  3.706693351 0.41924775  2.8849828582
#> 6                      fyear1994 -0.118952202 0.40958418 -0.9217224479
#> 7                      fyear1995  0.043060162 0.40708786 -0.7548173894
#> 8                      fyear1996  0.100803326 0.41182762 -0.7063639814
#> 9                      fyear1997 -0.595877590 0.41349979 -1.4063222959
#> 10                     fyear1998 -0.731133573 0.40940483 -1.5335522945
#> 11                     fyear1999 -0.418026699 0.40426484 -1.2103712192
#> 12                     fyear2000 -0.008214615 0.42492808 -0.8410583542
#> 13                     fyear2001  0.137378708 0.40808148 -0.6624463029
#> 14                     fyear2002  0.408813787 0.41429670 -0.4031928310
#> 15                     fyear2003  0.246142842 0.41430827 -0.5658864478
#> 16                     fyear2004 -0.969159441 0.40269880 -1.7584345819
#> 17                     fyear2005  0.025554985 0.39186036 -0.7424771997
#> 18                     fyear2006  0.602214237 0.39269289 -0.1674496787
#> 19                     fyear2007  0.351034556 0.38921615 -0.4118150731
#> 20                     fyear2008  0.301558233 0.39036135 -0.4635359630
#> 21                     fyear2009  0.080400536 0.38958103 -0.6831642497
#> 22                     fyear2010  0.662656166 0.39188814 -0.1054304685
#> 23                     fyear2011  0.247734357 0.39512575 -0.5266978816
#> 24                     fyear2012  0.100854403 0.39259612 -0.6686198511
#> 25                     fyear2013  0.466108433 0.38842681 -0.2951941212
#> 26                     fyear2014  0.156389468 0.39370883 -0.6152656656
#> 27                     fyear2015 -0.267006737 0.39316519 -1.0375963409
#> 28                     fyear2016 -0.151282214 0.39280337 -0.9211626677
#> 29                     fyear2017 -0.115556315 0.39224170 -0.8843359205
#> 30                     fyear2018  0.238284238 0.39501972 -0.5359401781
#> 31                     fyear2019 -0.186458600 0.42712577 -1.0236097344
#> 32                     fyear2020 -0.327708208 0.41651490 -1.1440624056
#> 33                      depth_sc  0.097959908 0.04974479  0.0004619065
#> 34                   depth_sq_sc -0.055185737 0.02318019 -0.1006180743
#> 35                       temp_sc  0.837214787 0.08491743  0.6707796801
#> 36                    temp_sq_sc -0.009634014 0.04451165 -0.0968752532
#> 37                        oxy_sc  0.563365911 0.05032856  0.4647237384
#> 38                        sal_sc  0.357375668 0.07569706  0.2090121600
#>       conf.high
#> 1   5.214642638
#> 2   4.613397247
#> 3   4.546060712
#> 4   4.744145232
#> 5   4.528403844
#> 6   0.683818045
#> 7   0.840937714
#> 8   0.907970633
#> 9   0.214567115
#> 10  0.071285149
#> 11  0.374317821
#> 12  0.824629125
#> 13  0.937203718
#> 14  1.220820404
#> 15  1.058172133
#> 16 -0.179884301
#> 17  0.793587169
#> 18  1.371878154
#> 19  1.113884186
#> 20  1.066652430
#> 21  0.843965321
#> 22  1.430742801
#> 23  1.022166596
#> 24  0.870328657
#> 25  1.227410987
#> 26  0.928044602
#> 27  0.503582867
#> 28  0.618598239
#> 29  0.653223290
#> 30  1.012508654
#> 31  0.650692535
#> 32  0.488645990
#> 33  0.195457910
#> 34 -0.009753399
#> 35  1.003649893
#> 36  0.077607226
#> 37  0.662008083
#> 38  0.505739176

q4

# Small cod model q4 spatial
mcod_s_q4_space <- sdmTMB(small_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                            temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                          data = filter(cpue2, quarter == 4),
                          mesh = spde_q4,
                          family = tweedie(link = "log"),
                          spatiotemporal = "IID",
                          spatial = "on",
                          time = "year",
                          reml = FALSE,
                          control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 5,271 rows (59%), 3,595 rows remaining

sanity(mcod_s_q4_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_s_q4_space, conf.int = TRUE)
#>                             term    estimate  std.error    conf.low   conf.high
#> 1              fsubstratebedrock  0.09435785 0.74557893 -1.36694999  1.55566570
#> 2            fsubstratehard clay -0.78082520 0.55629295 -1.87113935  0.30948894
#> 3  fsubstratehard-bottom complex -1.30971300 0.57320481 -2.43317377 -0.18625222
#> 4                  fsubstratemud -1.01679675 0.55400257 -2.10262183  0.06902833
#> 5                 fsubstratesand -0.81740941 0.55604348 -1.90723460  0.27241578
#> 6                      fyear1994  0.09575997 0.69520069 -1.26680834  1.45832828
#> 7                      fyear1995  1.08520169 0.65733901 -0.20315911  2.37356248
#> 8                      fyear1996 -0.65602062 0.69879535 -2.02563434  0.71359309
#> 9                      fyear1997  0.70967734 0.63506435 -0.53502593  1.95438060
#> 10                     fyear1998  0.42678207 0.65481121 -0.85662432  1.71018847
#> 11                     fyear1999  0.70832897 0.63583473 -0.53788420  1.95454215
#> 12                     fyear2000  1.15887477 0.61295487 -0.04249471  2.36024425
#> 13                     fyear2001  1.94399364 0.60273914  0.76264662  3.12534065
#> 14                     fyear2002  1.48091157 0.60100437  0.30296466  2.65885849
#> 15                     fyear2003 -0.17510629 0.61432459 -1.37916036  1.02894778
#> 16                     fyear2004  2.43321357 0.59484838  1.26733217  3.59909496
#> 17                     fyear2005  1.51499545 0.57814968  0.38184290  2.64814800
#> 18                     fyear2006  1.91959803 0.57576174  0.79112576  3.04807031
#> 19                     fyear2007  2.22121501 0.57216227  1.09979756  3.34263245
#> 20                     fyear2008  1.67400570 0.57090638  0.55504976  2.79296165
#> 21                     fyear2009  2.37379404 0.56946394  1.25766523  3.48992284
#> 22                     fyear2010  1.72512674 0.57286923  0.60232369  2.84792979
#> 23                     fyear2011  0.87635038 0.57624420 -0.25306751  2.00576826
#> 24                     fyear2012  1.31752499 0.58228856  0.17626038  2.45878960
#> 25                     fyear2013  1.98309676 0.57388746  0.85829801  3.10789550
#> 26                     fyear2014  0.88658990 0.58050683 -0.25118259  2.02436238
#> 27                     fyear2015  0.50716436 0.58397545 -0.63740649  1.65173521
#> 28                     fyear2016  0.16799697 0.58158813 -0.97189482  1.30788876
#> 29                     fyear2017  1.49569425 0.57485496  0.36899923  2.62238927
#> 30                     fyear2018  0.53595781 0.57927770 -0.59940563  1.67132124
#> 31                     fyear2019  0.05818354 0.64153030 -1.19919274  1.31555983
#> 32                     fyear2020  0.26687992 0.61174312 -0.93211457  1.46587441
#> 33                      depth_sc  0.17234560 0.09341009 -0.01073482  0.35542602
#> 34                   depth_sq_sc -1.27936173 0.08813100 -1.45209531 -1.10662815
#> 35                       temp_sc  0.76108938 0.13370443  0.49903352  1.02314524
#> 36                    temp_sq_sc -0.20875558 0.07211844 -0.35010512 -0.06740604
#> 37                        oxy_sc  1.07125991 0.10102534  0.87325388  1.26926595
#> 38                        sal_sc  0.44360031 0.11241137  0.22327807  0.66392255

# Large cod model q4 spatial
mcod_l_q4_space <- sdmTMB(large_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                            temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                          data = filter(cpue2, quarter == 4),
                          mesh = spde_q4,
                          family = tweedie(link = "log"),
                          spatiotemporal = "IID",
                          spatial = "on",
                          time = "year",
                          reml = FALSE,
                          control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 5,271 rows (59%), 3,595 rows remaining

sanity(mcod_l_q4_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_l_q4_space, conf.int = TRUE)
#>                             term    estimate  std.error    conf.low   conf.high
#> 1              fsubstratebedrock  4.21170576 0.68969740  2.85992369  5.56348782
#> 2            fsubstratehard clay  3.86010397 0.51098291  2.85859587  4.86161206
#> 3  fsubstratehard-bottom complex  3.38036271 0.52245906  2.35636178  4.40436364
#> 4                  fsubstratemud  3.77126497 0.50969618  2.77227881  4.77025113
#> 5                 fsubstratesand  3.71416028 0.51148436  2.71166936  4.71665120
#> 6                      fyear1994  0.42759184 0.56627641 -0.68228954  1.53747321
#> 7                      fyear1995  1.17104458 0.56088597  0.07172829  2.27036088
#> 8                      fyear1996  0.30666879 0.57086084 -0.81219790  1.42553548
#> 9                      fyear1997 -0.27644075 0.55060815 -1.35561290  0.80273140
#> 10                     fyear1998 -0.12468170 0.55891250 -1.22013007  0.97076667
#> 11                     fyear1999 -0.50545258 0.55787366 -1.59886486  0.58795970
#> 12                     fyear2000  0.16521615 0.53777897 -0.88881126  1.21924356
#> 13                     fyear2001  0.17695597 0.53395591 -0.86957838  1.22349032
#> 14                     fyear2002  0.82122342 0.52673599 -0.21116016  1.85360699
#> 15                     fyear2003 -1.41070205 0.54015717 -2.46939066 -0.35201344
#> 16                     fyear2004  0.11698731 0.53032500 -0.92243058  1.15640520
#> 17                     fyear2005  0.60911272 0.51052745 -0.39150268  1.60972813
#> 18                     fyear2006  0.42784518 0.51125268 -0.57419165  1.42988202
#> 19                     fyear2007  0.57901828 0.50949887 -0.41958116  1.57761772
#> 20                     fyear2008  0.71078769 0.50633284 -0.28160644  1.70318182
#> 21                     fyear2009  0.48858518 0.50727701 -0.50565948  1.48282984
#> 22                     fyear2010  0.74756835 0.50853799 -0.24914779  1.74428449
#> 23                     fyear2011  0.01460874 0.51048187 -0.98591734  1.01513482
#> 24                     fyear2012 -0.40274524 0.51964102 -1.42122292  0.61573244
#> 25                     fyear2013  0.06254185 0.51097323 -0.93894727  1.06403097
#> 26                     fyear2014 -0.13121984 0.51357630 -1.13781090  0.87537122
#> 27                     fyear2015 -0.02210225 0.51330778 -1.02816700  0.98396251
#> 28                     fyear2016 -0.71613915 0.51517454 -1.72586269  0.29358439
#> 29                     fyear2017 -0.41320127 0.51291329 -1.41849284  0.59209030
#> 30                     fyear2018 -1.08808137 0.51586811 -2.09916429 -0.07699845
#> 31                     fyear2019 -1.21055219 0.56145868 -2.31099099 -0.11011339
#> 32                     fyear2020 -1.50211588 0.54522694 -2.57074105 -0.43349071
#> 33                      depth_sc  0.37248818 0.07157203  0.23220958  0.51276677
#> 34                   depth_sq_sc -0.57573368 0.05693097 -0.68731633 -0.46415102
#> 35                       temp_sc  0.28283821 0.09928947  0.08823443  0.47744200
#> 36                    temp_sq_sc -0.08858374 0.05410450 -0.19462660  0.01745913
#> 37                        oxy_sc  0.73749983 0.07553416  0.58945561  0.88554405
#> 38                        sal_sc  0.24195972 0.08730570  0.07084370  0.41307574

# Flounder model q4 spatial
mfle_q4_space <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                          temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                      data = filter(cpue2, quarter == 4),
                      mesh = spde_q4,
                      family = tweedie(link = "log"),
                      spatiotemporal = "IID",
                      spatial = "on",
                      time = "year",
                      reml = FALSE,
                      control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 5,271 rows (59%), 3,595 rows remaining

sanity(mfle_q4_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfle_q4_space, conf.int = TRUE)
#>                             term    estimate  std.error    conf.low   conf.high
#> 1              fsubstratebedrock  2.51946863 0.55705485  1.42766119  3.61127606
#> 2            fsubstratehard clay  2.15626423 0.51660542  1.14373621  3.16879225
#> 3  fsubstratehard-bottom complex  2.00553041 0.52214542  0.98214418  3.02891664
#> 4                  fsubstratemud  2.09766143 0.51521436  1.08785984  3.10746303
#> 5                 fsubstratesand  2.17656896 0.51672767  1.16380134  3.18933657
#> 6                      fyear1994 -0.47801146 0.56314064 -1.58174682  0.62572391
#> 7                      fyear1995  0.70747231 0.54561884 -0.36192096  1.77686558
#> 8                      fyear1996  0.21787202 0.56110562 -0.88187479  1.31761884
#> 9                      fyear1997  0.21476366 0.52933021 -0.82270450  1.25223181
#> 10                     fyear1998 -0.27099120 0.53917133 -1.32774758  0.78576518
#> 11                     fyear1999 -1.03874014 0.54195660 -2.10095556  0.02347528
#> 12                     fyear2000  0.02358672 0.52036980 -0.99631935  1.04349279
#> 13                     fyear2001 -0.13931059 0.51746380 -1.15352100  0.87489981
#> 14                     fyear2002  0.79548001 0.50897871 -0.20209993  1.79305996
#> 15                     fyear2003 -1.03066836 0.50997714 -2.03020519 -0.03113152
#> 16                     fyear2004  0.39486395 0.50887822 -0.60251904  1.39224694
#> 17                     fyear2005  0.43871414 0.49061216 -0.52286802  1.40029631
#> 18                     fyear2006  0.46242684 0.48906494 -0.49612283  1.42097650
#> 19                     fyear2007  0.59181644 0.48775073 -0.36415742  1.54779031
#> 20                     fyear2008  0.80776983 0.48502260 -0.14285700  1.75839667
#> 21                     fyear2009  0.71303600 0.48347509 -0.23455776  1.66062976
#> 22                     fyear2010  1.20423026 0.48405021  0.25550928  2.15295124
#> 23                     fyear2011  0.75829681 0.48413917 -0.19059854  1.70719215
#> 24                     fyear2012  0.36522721 0.48772824 -0.59070257  1.32115699
#> 25                     fyear2013  0.68792578 0.48757151 -0.26769682  1.64354837
#> 26                     fyear2014  0.34406155 0.48797468 -0.61235125  1.30047435
#> 27                     fyear2015  0.37608552 0.48740447 -0.57920969  1.33138074
#> 28                     fyear2016  0.25798676 0.48543592 -0.69345017  1.20942369
#> 29                     fyear2017  0.47267936 0.48773506 -0.48326379  1.42862251
#> 30                     fyear2018  0.25476915 0.48878170 -0.70322538  1.21276368
#> 31                     fyear2019  0.26560882 0.52792617 -0.76910747  1.30032511
#> 32                     fyear2020 -0.12599992 0.51004667 -1.12567302  0.87367319
#> 33                      depth_sc -0.18628399 0.06551021 -0.31468164 -0.05788635
#> 34                   depth_sq_sc -0.70612932 0.06538288 -0.83427741 -0.57798124
#> 35                       temp_sc  0.64874416 0.09993836  0.45286857  0.84461975
#> 36                    temp_sq_sc -0.06611203 0.05554101 -0.17497041  0.04274636
#> 37                        oxy_sc  1.20548898 0.07915864  1.05034089  1.36063707
#> 38                        sal_sc  0.22610329 0.08758623  0.05443743  0.39776914

Plot coefficients

# Q1
# small cod
small_cod_q1 <- tidy(mcod_s_q1_space, effects = "fixed", conf.int = TRUE) %>%
  mutate(species = "Small cod",
         quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# large cod
large_cod_q1 <- tidy(mcod_l_q1_space, effects = "fixed", conf.int = TRUE) %>%
  mutate(species = "Large cod",
         quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# flounder
flounder_q1 <- tidy(mfle_q1_space, effects = "fixed", conf.int = TRUE) %>%
  mutate(species = "Flounder",
         quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# Q4
# small cod
small_cod_q4 <- tidy(mcod_s_q4_space, effects = "fixed", conf.int = TRUE) %>%
  mutate(species = "Small cod",
         quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# large cod
large_cod_q4 <- tidy(mcod_l_q4_space, effects = "fixed", conf.int = TRUE) %>%
  mutate(species = "Large cod",
         quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# flounder
flounder_q4 <- tidy(mfle_q4_space, effects = "fixed", conf.int = TRUE) %>%
  mutate(species = "Flounder",
         quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA


coefs <- bind_rows(small_cod_q1,
                   large_cod_q1,
                   flounder_q1,
                      
                   small_cod_q4,
                   large_cod_q4,
                   flounder_q4) %>% 
  mutate(param_group = ifelse(term %in% c("fsubstratebedrock",
                                          "fsubstratehard clay",
                                          "fsubstratehard-bottom complex",
                                          "fsubstratemud",
                                          "fsubstratesand"),
                              "Substrate",
                              "Continious")) %>% 
  mutate(term = ifelse(term == "fsubstratebedrock", "Bedrock", term),
         term = ifelse(term == "fsubstratehard clay", "Hard clay", term),
         term = ifelse(term == "fsubstratehard-bottom complex", "Hard-bottom\ncomplex", term),
         term = ifelse(term == "fsubstratemud", "Mud", term),
         term = ifelse(term == "fsubstratesand", "Sand", term),
         term = ifelse(term == "depth_sc", "Depth", term),
         term = ifelse(term == "depth_sq_sc", "Depth^2", term),
         term = ifelse(term == "temp_sc", "Temperature", term),
         term = ifelse(term == "temp_sq_sc", "Temperature^2", term),
         term = ifelse(term == "oxy_sc", "Oxygen", term),
         term = ifelse(term == "sal_sc", "Salinity", term))
#> mutate: new variable 'param_group' (character) with 2 unique values and 0% NA
#> mutate: changed 66 values (29%) of 'term' (0 new NA)

coefs %>% 
  filter(!grepl('year', term)) %>% 
  ggplot(aes(term, estimate, color = species, shape = factor(quarter))) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.2) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  labs(shape = "Quarter") + 
  facet_wrap(~param_group, scales = "free") +
  labs(x = "Predictor", y = "Standardized coefficient") +
  theme_plot() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90),
        aspect.ratio = 1) +
  NULL 
#> filter: removed 162 rows (71%), 66 rows remaining

ggsave("figures/habitat_coefs.pdf", width = 20, height = 20, units = "cm")

Predict on grid

mcod_s_q1_pred <- predict(mcod_s_q1_space, newdata = pred_grid_q1)
mcod_s_q4_pred <- predict(mcod_s_q4_space, newdata = pred_grid_q4)

mcod_l_q1_pred <- predict(mcod_l_q1_space, newdata = pred_grid_q1)
mcod_l_q4_pred <- predict(mcod_l_q4_space, newdata = pred_grid_q4)

mfle_q1_pred <- predict(mfle_q1_space, newdata = pred_grid_q1)
mfle_q4_pred <- predict(mfle_q4_space, newdata = pred_grid_q4)

Plot predictions on map

p1 <- plot_map + 
  geom_raster(data = bind_rows(mcod_s_q1_pred, mcod_s_q4_pred) %>%
                filter(year %in% c("1995", "2019")),
              aes(X*1000, Y*1000, fill = exp(est))) +
  facet_grid(quarter~year) +
  ggtitle("Cod < 25 cm") +
  scale_fill_viridis(trans = "sqrt")
#> filter: removed 550,420 rows (93%), 42,340 rows remaining

p2 <- plot_map + 
  geom_raster(data = bind_rows(mcod_l_q1_pred, mcod_l_q4_pred) %>%
                filter(year %in% c("1995", "2019")),
              aes(X*1000, Y*1000, fill = exp(est))) +
  facet_grid(quarter~year) +
  ggtitle("Cod > 25 cm") +
  scale_fill_viridis(trans = "sqrt")
#> filter: removed 550,420 rows (93%), 42,340 rows remaining

p3 <- plot_map + 
  geom_raster(data = bind_rows(mfle_q1_pred, mfle_q4_pred) %>%
                filter(year %in% c("1995", "2019")),
              aes(X*1000, Y*1000, fill = exp(est))) +
  facet_grid(quarter~year) +
  ggtitle("Flounder") +
  scale_fill_viridis(trans = "sqrt")
#> filter: removed 550,420 rows (93%), 42,340 rows remaining

((p1 + p2) / (p3 + plot_spacer())) + plot_annotation(tag_levels = "A") &
  theme(legend.text = element_text(angle = 90))

# Check flounder...

plot_map + 
  geom_raster(data = mfle_q1_pred,
              aes(X*1000, Y*1000, fill = exp(est))) +
  facet_wrap(~year) +
  scale_fill_viridis(trans = "sqrt")

plot_map + 
  geom_raster(data = mfle_q4_pred,
              aes(X*1000, Y*1000, fill = exp(est))) +
  facet_wrap(~year) +
  scale_fill_viridis(trans = "sqrt")

Calculate spatial overlap indices

all_pred_q1 <- mcod_s_q1_pred %>% mutate(small_cod = exp(est))
#> mutate: new variable 'small_cod' (double) with 296,380 unique values and 0% NA

all_pred_q1 <- all_pred_q1 %>%
  mutate(large_cod = exp(mcod_l_q1_pred$est),
         fle = exp(mfle_q1_pred$est))
#> mutate: new variable 'large_cod' (double) with 296,380 unique values and 0% NA
#>         new variable 'fle' (double) with 296,380 unique values and 0% NA

all_pred_q4 <- mcod_s_q4_pred %>% mutate(small_cod = exp(est))
#> mutate: new variable 'small_cod' (double) with 296,380 unique values and 0% NA

all_pred_q4 <- all_pred_q4 %>%
  mutate(large_cod = exp(mcod_l_q4_pred$est),
         fle = exp(mfle_q4_pred$est))
#> mutate: new variable 'large_cod' (double) with 296,380 unique values and 0% NA
#>         new variable 'fle' (double) with 296,380 unique values and 0% NA

all_pred <- bind_rows(all_pred_q1, all_pred_q4)

# Calculate biomass overlap
loc_collocspfn <- function(prey, pred) {
  p_prey <- prey/sum(prey, na.rm = T)
  p_pred <- pred/sum(pred, na.rm = T)
  (p_prey*p_pred)/(sqrt(sum(p_prey^2, na.rm = T))* sqrt(sum(p_pred^2, na.rm = T)))
}

loc_collocspfn_tot <- function(prey, pred) {
  p_prey <- prey/sum(prey, na.rm = T)
  p_pred <- pred/sum(pred, na.rm = T)
  sum((p_prey*p_pred))/(sqrt(sum(p_prey^2, na.rm = T))* sqrt(sum(p_pred^2, na.rm = T)))
}

all_pred <- all_pred %>% 
  group_by(year, quarter) %>% 
  mutate(scod_fle_ovr = loc_collocspfn(pred = small_cod, prey = fle),
         scod_fle_ovr_tot = loc_collocspfn_tot(pred = small_cod, prey = fle),
         lcod_fle_ovr = loc_collocspfn(pred = large_cod, prey = fle),
         lcod_fle_ovr_tot = loc_collocspfn_tot(pred = large_cod, prey = fle)) %>% 
  ungroup()
#> group_by: 2 grouping variables (year, quarter)
#> mutate (grouped): new variable 'scod_fle_ovr' (double) with 592,760 unique values and 0% NA
#>                   new variable 'scod_fle_ovr_tot' (double) with 56 unique values and 0% NA
#>                   new variable 'lcod_fle_ovr' (double) with 592,760 unique values and 0% NA
#>                   new variable 'lcod_fle_ovr_tot' (double) with 56 unique values and 0% NA
#> ungroup: no grouping variables

Plot

plot_map_fc +
  geom_raster(data = all_pred %>% filter(quarter == 1),
              aes(x = X, y = Y, fill = scod_fle_ovr)) +
  scale_fill_viridis_c(trans = "sqrt", name = "small cod-flounder") +
  facet_wrap(~ year, ncol = 5) +
  theme(legend.text = element_text(angle = 90)) + 
  NULL
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
#> Warning: Removed 296380 rows containing missing values (geom_raster).

plot_map_fc +
  geom_raster(data = all_pred %>% filter(quarter == 4),
              aes(x = X, y = Y, fill = lcod_fle_ovr)) +
  scale_fill_viridis_c(trans = "sqrt", name = "large cod-flounder") +
  facet_wrap(~ year, ncol = 5) +
  theme(legend.text = element_text(angle = 90)) + 
  NULL
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
#> Warning: Removed 296380 rows containing missing values (geom_raster).

# In space for selected year
all_pred2 <- all_pred %>%
  filter(year %in% c(1995, 2015)) %>% 
  dplyr::select(X, Y, scod_fle_ovr, lcod_fle_ovr, quarter, year) %>% 
  pivot_longer(c(scod_fle_ovr, lcod_fle_ovr)) %>% 
  rename(overlap = name)
#> filter: removed 550,420 rows (93%), 42,340 rows remaining
#> pivot_longer: reorganized (scod_fle_ovr, lcod_fle_ovr) into (name, value) [was 42340x6, now 84680x6]
#> rename: renamed one variable (overlap)

p1 <- plot_map_fc +
  geom_raster(data = filter(all_pred2, overlap == "scod_fle_ovr"),
              aes(x = X*1000, y = Y*1000, fill = value)) +
  scale_fill_viridis_c(trans = "sqrt", name = "Overlap") +
  facet_grid(year ~ quarter) +
  theme(legend.text = element_text(angle = 90)) + 
  ggtitle("Small cod : Flonder") +
  NULL
#> filter: removed 42,340 rows (50%), 42,340 rows remaining

p2 <- plot_map_fc +
  geom_raster(data = filter(all_pred2, overlap == "lcod_fle_ovr"),
              aes(x = X*1000, y = Y*1000, fill = value)) +
  scale_fill_viridis_c(trans = "sqrt", name = "Overlap") +
  facet_grid(year ~ quarter) +
  theme(legend.text = element_text(angle = 90)) + 
  ggtitle("Large cod : Flonder") +
  NULL
#> filter: removed 42,340 rows (50%), 42,340 rows remaining

p1 + p2

ggsave("figures/spatial_overlap.pdf", width = 22, height = 14, units = "cm")

# Over time
p5 <- all_pred %>% 
  distinct(year, quarter, .keep_all = TRUE) %>% 
  ggplot(aes(year, lcod_fle_ovr_tot, color = factor(quarter), fill = factor(quarter))) + 
  labs(color = "Quarter", x = "Year", y = "Overlap index") + 
  guides(fill = "none") +
  stat_smooth(method = "gam", formula = y~s(x, k = 3)) +
  ggtitle("Large cod : Flounder") +
  geom_point()
#> distinct: removed 592,704 rows (>99%), 56 rows remaining

p6 <- all_pred %>% 
  distinct(year, quarter, .keep_all = TRUE) %>% 
  ggplot(aes(year, scod_fle_ovr_tot, color = factor(quarter), fill = factor(quarter))) + 
  labs(color = "Quarter", x = "Year", y = "Overlap index") + 
  guides(fill = "none") +
  stat_smooth(method = "gam", formula = y~s(x, k = 3)) +
  ggtitle("Small cod : Flounder") +
  geom_point()
#> distinct: removed 592,704 rows (>99%), 56 rows remaining

(p5 | p6) + plot_annotation(tag_levels = "A") & theme(aspect.ratio = 1)

ggsave("figures/spatial_overlap_time.pdf", width = 20, height = 20, units = "cm")

Conditional effects

# Scale with respect to data!
nd_depth <- data.frame(depth = seq(min(cpue2$depth), max(cpue2$depth), length.out = 100)) %>% 
  mutate(depth_ct = depth - mean(cpue2$depth),
         depth_sq = depth_ct * depth_ct,
         depth_sq_sc = (depth_sq - mean(cpue2$depth_sq)) / sd(cpue2$depth_sq),
         depth_sc = (depth - mean(cpue2$depth)) / sd(cpue2$depth),
         temp_sq_sc = 0,
         temp_sc = 0,
         oxy_sc = 0,
         sal_sc = 0,
         year = 1999L,
         fyear = as.factor(1999),
         fsubstrate = "mud")
#> mutate: new variable 'depth_ct' (double) with 100 unique values and 0% NA
#>         new variable 'depth_sq' (double) with 100 unique values and 0% NA
#>         new variable 'depth_sq_sc' (double) with 100 unique values and 0% NA
#>         new variable 'depth_sc' (double) with 100 unique values and 0% NA
#>         new variable 'temp_sq_sc' (double) with one unique value and 0% NA
#>         new variable 'temp_sc' (double) with one unique value and 0% NA
#>         new variable 'oxy_sc' (double) with one unique value and 0% NA
#>         new variable 'sal_sc' (double) with one unique value and 0% NA
#>         new variable 'year' (integer) with one unique value and 0% NA
#>         new variable 'fyear' (factor) with one unique value and 0% NA
#>         new variable 'fsubstrate' (character) with one unique value and 0% NA

# Q1
margin_small_cod_pred_q1_depth <- predict(mcod_s_q1_space,
                                          newdata = nd_depth,
                                          se_fit = TRUE,
                                          re_form = NA) %>%
  mutate(species = "small_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_large_cod_pred_q1_depth <- predict(mcod_l_q1_space,
                                          newdata = nd_depth,
                                          se_fit = TRUE,
                                          re_form = NA) %>%
  mutate(species ="large_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_fle_pred_q1_depth <- predict(mfle_q1_space,
                                    newdata = nd_depth,
                                    se_fit = TRUE,
                                    re_form = NA) %>%
  mutate(species = "flounder", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# Q4
margin_small_cod_pred_q4_depth <- predict(mcod_s_q4_space,
                                          newdata = nd_depth,
                                          se_fit = TRUE,
                                          re_form = NA) %>%
  mutate(species = "small_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_large_cod_pred_q4_depth <- predict(mcod_l_q4_space,
                                          newdata = nd_depth,
                                          se_fit = TRUE,
                                          re_form = NA) %>%
  mutate(species = "large_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_fle_pred_q4_depth <- predict(mfle_q4_space,
                                    newdata = nd_depth,
                                    se_fit = TRUE,
                                    re_form = NA) %>%
  mutate(species = "flounder", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margins_depth <- bind_rows(margin_small_cod_pred_q1_depth,
                           margin_large_cod_pred_q1_depth,
                           margin_fle_pred_q1_depth,
                           margin_small_cod_pred_q4_depth,
                           margin_large_cod_pred_q4_depth,
                           margin_fle_pred_q4_depth) %>% 
    mutate(species = ifelse(species == "flounder", "Flounder", species),
           species = ifelse(species == "small_cod", "Small cod", species),
           species = ifelse(species == "large_cod", "Large cod", species))
#> mutate: changed 600 values (100%) of 'species' (0 new NA)

# Temperature
nd_temp <- data.frame(temp = seq(min(cpue2$temp), max(cpue2$temp), length.out = 100)) %>% 
  mutate(temp_ct = temp - mean(cpue2$temp),
         temp_sq = temp_ct * temp_ct,
         temp_sq_sc = (temp_sq - mean(cpue2$temp_sq)) / sd(cpue2$temp_sq),
         temp_sc = (temp - mean(cpue2$temp)) / sd(cpue2$temp),
         depth_sq_sc = 0,
         depth_sc = 0,
         oxy_sc = 0,
         sal_sc = 0,
         year = 1999L,
         fyear = as.factor(1999),
         fsubstrate = "mud")
#> mutate: new variable 'temp_ct' (double) with 100 unique values and 0% NA
#>         new variable 'temp_sq' (double) with 100 unique values and 0% NA
#>         new variable 'temp_sq_sc' (double) with 100 unique values and 0% NA
#>         new variable 'temp_sc' (double) with 100 unique values and 0% NA
#>         new variable 'depth_sq_sc' (double) with one unique value and 0% NA
#>         new variable 'depth_sc' (double) with one unique value and 0% NA
#>         new variable 'oxy_sc' (double) with one unique value and 0% NA
#>         new variable 'sal_sc' (double) with one unique value and 0% NA
#>         new variable 'year' (integer) with one unique value and 0% NA
#>         new variable 'fyear' (factor) with one unique value and 0% NA
#>         new variable 'fsubstrate' (character) with one unique value and 0% NA

# Q1
margin_small_cod_pred_q1_temp <- predict(mcod_s_q1_space,
                                         newdata = nd_temp,
                                         se_fit = TRUE,
                                         re_form = NA) %>%
  mutate(species = "small_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_large_cod_pred_q1_temp <- predict(mcod_l_q1_space,
                                         newdata = nd_temp,
                                         se_fit = TRUE,
                                         re_form = NA) %>%
  mutate(species ="large_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_fle_pred_q1_temp <- predict(mfle_q1_space,
                                   newdata = nd_temp,
                                   se_fit = TRUE,
                                   re_form = NA) %>%
  mutate(species = "flounder", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# Q4
margin_small_cod_pred_q4_temp <- predict(mcod_s_q4_space,
                                         newdata = nd_temp,
                                         se_fit = TRUE,
                                         re_form = NA) %>%
  mutate(species = "small_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_large_cod_pred_q4_temp <- predict(mcod_l_q4_space,
                                         newdata = nd_temp,
                                         se_fit = TRUE,
                                         re_form = NA) %>%
  mutate(species = "large_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_fle_pred_q4_temp <- predict(mfle_q4_space,
                                   newdata = nd_temp,
                                   se_fit = TRUE,
                                   re_form = NA) %>%
  mutate(species = "flounder", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margins_temp <- bind_rows(margin_small_cod_pred_q1_temp,
                          margin_large_cod_pred_q1_temp,
                          margin_fle_pred_q1_temp,
                          margin_small_cod_pred_q4_temp,
                          margin_large_cod_pred_q4_temp,
                          margin_fle_pred_q4_temp) %>% 
    mutate(species = ifelse(species == "flounder", "Flounder", species),
           species = ifelse(species == "small_cod", "Small cod", species),
           species = ifelse(species == "large_cod", "Large cod", species))
#> mutate: changed 600 values (100%) of 'species' (0 new NA)

# Plot!
p1 <- ggplot(margins_depth, 
             aes(depth, exp(est), ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se),
                 color = species, fill = species)) +
  geom_line() +
  facet_wrap(~quarter, scales = "free") +
  geom_ribbon(alpha = 0.4, color = NA) +
  coord_cartesian(xlim = c(10, 170)) +
  theme(aspect.ratio = 1) +
  labs(x = "Depth (m)", y = "Biomass density", color = "", fill = "") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)


p2 <- ggplot(margins_temp,
             aes(temp, exp(est), ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se),
                 color = species, fill = species)) +
  geom_line() +
  facet_wrap(~quarter, scales = "free") +
  geom_ribbon(alpha = 0.4, color = NA) +
  coord_cartesian(xlim = c(1, 14)) +
  theme(aspect.ratio = 1) +
  labs(x = "Temperature (°C)", y = "Biomass density", color = "", fill = "") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)

(p1 / p2) + plot_annotation(tag_levels = "A") + plot_layout(guides = "collect")

ggsave("figures/supp/conditional_effects.pdf", width = 20, height = 20, units = "cm")

# Scaled
p3 <- margins_temp %>% 
  ungroup() %>% 
  group_by(species, quarter) %>% 
  mutate(max_est = max(exp(est)),
         est_sc = exp(est) / max_est) %>% 
  ggplot(aes(temp, est_sc, color = species, fill = species)) +
  geom_line(size = 1.3) +
  facet_wrap(~quarter, scales = "free") +
  coord_cartesian(xlim = c(1, 14)) +
  theme(aspect.ratio = 1) +
  labs(x = "Temperature (°C)", y = "Scaled biomass density", color = "", fill = "") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#> ungroup: no grouping variables
#> group_by: 2 grouping variables (species, quarter)
#> mutate (grouped): new variable 'max_est' (double) with 6 unique values and 0% NA
#>                   new variable 'est_sc' (double) with 595 unique values and 0% NA

p4 <- margins_depth %>% 
  ungroup() %>% 
  group_by(species, quarter) %>% 
  mutate(max_est = max(exp(est)),
         est_sc = exp(est) / max_est) %>% 
  ggplot(aes(depth, est_sc, color = species, fill = species)) +
  geom_line(size = 1.3) +
  facet_wrap(~quarter, scales = "free") +
  coord_cartesian(xlim = c(10, 170)) +
  theme(aspect.ratio = 1) +
  labs(x = "Depth (m)", y = "Scaled biomass density", color = "", fill = "") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#> ungroup: no grouping variables
#> group_by: 2 grouping variables (species, quarter)
#> mutate (grouped): new variable 'max_est' (double) with 6 unique values and 0% NA
#>                   new variable 'est_sc' (double) with 595 unique values and 0% NA

(p3 / p4) + plot_annotation(tag_levels = "A") + plot_layout(guides = "collect")

ggsave("figures/conditional_effects.pdf", width = 20, height = 20, units = "cm")

Weighted quantiles

head(all_pred) %>% as.data.frame()
#>     X    Y depth year      lon      lat substrate quarter      oxy     temp
#> 1 450 5984    11 1993 14.23718 54.00188      sand       1 8.713132 2.821855
#> 2 454 5984    11 1993 14.29820 54.00225      sand       1 8.743857 2.811376
#> 3 458 5984    11 1993 14.35922 54.00259      sand       1 8.745697 2.806495
#> 4 462 5984    12 1993 14.42024 54.00290      sand       1 8.759304 2.802808
#> 5 466 5984    12 1993 14.48127 54.00318      sand       1 8.765212 2.805644
#> 6 470 5984    11 1993 14.54229 54.00343      sand       1 8.805174 2.792613
#>        sal sub_div ices_rect density_saduria biomass_spr biomass_her
#> 1 7.293935      24      37G4               0          NA          NA
#> 2 7.346682      24      37G4               0          NA          NA
#> 3 7.383028      24      37G4               0          NA          NA
#> 4 7.418652      24      37G4               0          NA          NA
#> 5 7.501428      24      37G4               0          NA          NA
#> 6 7.585449      24      37G4               0          NA          NA
#>   biomass_spr_sd biomass_her_sd  depth_ct depth_sq depth_sq_sc  depth_sc
#> 1             NA             NA -36.63035 1341.783   0.4698617 -1.343238
#> 2             NA             NA -36.63035 1341.783   0.4698617 -1.343238
#> 3             NA             NA -36.63035 1341.783   0.4698617 -1.343238
#> 4             NA             NA -35.63035 1269.522   0.4131041 -1.306567
#> 5             NA             NA -35.63035 1269.522   0.4131041 -1.306567
#> 6             NA             NA -36.63035 1341.783   0.4698617 -1.343238
#>     temp_ct  temp_sq temp_sq_sc   temp_sc   oxy_sc    sal_sc fyear fsubstrate
#> 1 -3.484890 12.14446  0.5884181 -1.305532 1.178372 -1.098960  1993       sand
#> 2 -3.495369 12.21760  0.5969916 -1.309458 1.190285 -1.082916  1993       sand
#> 3 -3.500250 12.25175  0.6009941 -1.311287 1.190999 -1.071861  1993       sand
#> 4 -3.503937 12.27757  0.6040213 -1.312668 1.196275 -1.061026  1993       sand
#> 5 -3.501101 12.25771  0.6016927 -1.311606 1.198566 -1.035849  1993       sand
#> 6 -3.514132 12.34912  0.6124078 -1.316487 1.214061 -1.010293  1993       sand
#>         est est_non_rf     est_rf    omega_s epsilon_st  small_cod large_cod
#> 1 -3.656425  -3.277790 -0.3786354 -0.5617690  0.1831337 0.02582466  5.541801
#> 2 -3.865310  -3.278880 -0.5864303 -0.7847163  0.1982860 0.02095642  5.108481
#> 3 -4.067168  -3.282109 -0.7850593 -0.9984999  0.2134406 0.01712582  4.721942
#> 4 -4.041875  -3.169858 -0.8720176 -1.1006398  0.2286222 0.01756450  4.789478
#> 5 -4.125475  -3.166499 -0.9589759 -1.2027798  0.2438039 0.01615582  4.634219
#> 6 -4.325476  -3.279542 -1.0459342 -1.3049197  0.2589855 0.01322725  4.266585
#>        fle scod_fle_ovr scod_fle_ovr_tot lcod_fle_ovr lcod_fle_ovr_tot
#> 1 14.79964 4.175460e-08        0.3650466 9.478820e-07        0.5484333
#> 2 14.55492 3.332310e-08        0.3650466 8.593177e-07        0.5484333
#> 3 14.27225 2.670313e-08        0.3650466 7.788702e-07        0.5484333
#> 4 14.71945 2.824528e-08        0.3650466 8.147639e-07        0.5484333
#> 5 15.16271 2.676235e-08        0.3650466 8.120923e-07        0.5484333
#> 6 15.55304 2.247519e-08        0.3650466 7.669161e-07        0.5484333

all_pred_sub <- all_pred %>%
  dplyr::select(small_cod, large_cod, fle, year, temp, oxy, sal, depth, quarter) %>% 
  pivot_longer(c(small_cod, large_cod, fle)) %>% 
  rename(Species = name, density = value) 
#> pivot_longer: reorganized (small_cod, large_cod, fle) into (name, value) [was 592760x9, now 1778280x8]
#> rename: renamed 2 variables (Species, density)

wm <- all_pred_sub %>%
  pivot_longer(c(temp, oxy, sal, depth)) %>% 
  mutate(name = ifelse(name == "temp", "Temperature (°C)", name),
         name = ifelse(name == "sal", "Salinity", name),
         name = ifelse(name == "depth", "Depth", name),
         name = ifelse(name == "oxy", "Oxygen (ml/L)", name)) %>% 
  group_by(year, quarter, Species, name) %>%
  summarise("Weighted median" = hutils::weighted_quantile(v = value, w = density, p = c(0.5)),
            "1st quartile" = hutils::weighted_quantile(v = value, w = density, p = c(0.25)),
            "3rd quartile" = hutils::weighted_quantile(v = value, w = density, p = c(0.75))) %>% 
  rename(env_var = name) %>% 
  mutate(Species = ifelse(Species == "fle", "Flounder", Species),
         Species = ifelse(Species == "large_cod", "Large cod", Species),
         Species = ifelse(Species == "small_cod", "Small cod", Species))
#> pivot_longer: reorganized (temp, oxy, sal, depth) into (name, value) [was 1778280x8, now 7113120x6]
#> mutate: changed 7,113,120 values (100%) of 'name' (0 new NA)
#> group_by: 4 grouping variables (year, quarter, Species, name)
#> summarise: now 672 rows and 7 columns, 3 group variables remaining (year, quarter, Species)
#> rename: renamed one variable (env_var)
#> mutate (grouped): changed 672 values (100%) of 'Species' (0 new NA)

ggplot(wm, aes(year, `Weighted median`, color = Species, fill = Species)) +
  geom_line() + 
  ggh4x::facet_grid2(env_var ~ quarter, scales = "free_y", independent = "y") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)

ggplot(wm, aes(year, `Weighted median`, color = Species)) +
  geom_line(size = 0.8) + 
  labs(linetype = "Quarter", x = "Year") +
  ggh4x::facet_grid2(env_var ~ quarter, scales = "free_y", independent = "y") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)

ggsave("figures/weighted_quantiles.pdf", width = 20, height = 20, units = "cm")

# Scaled
wm %>% 
  group_by(Species, quarter, env_var) %>% 
  mutate(`Weighted median scaled` = `Weighted median` - mean(`Weighted median`)) %>% 
  ggplot(aes(year, `Weighted median scaled`, color = Species)) +
  geom_line(size = 0.8) + 
  labs(linetype = "Quarter", x = "Year") +
  ggh4x::facet_grid2(env_var ~ quarter, scales = "free_y", independent = "y") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#> group_by: 3 grouping variables (Species, quarter, env_var)
#> mutate (grouped): new variable 'Weighted median scaled' (double) with 647 unique values and 0% NA

ggsave("figures/supp/weighted_quantiles_scaled.pdf", width = 20, height = 20, units = "cm")

# ggplot(wm, aes(year, Median, color = Species, fill = Species)) +
#   geom_ribbon(aes(ymin = `1st quartile`, ymax = `3rd quartile`), alpha = 0.4, color = NA) + 
#   facet_wrap(env_var ~ quarter, scales = "free", ncol = 2) + 
#   scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
#   scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
# 
# ggsave("figures/supp/weighted_quantiles_range.pdf", width = 20, height = 20, units = "cm")

[Extra for bentfish - calculate biomass index per subdivision]

Compare with Q4 from previous index

cpue3 <- cpue2 %>%
  mutate(tot_cod = large_cod + small_cod) %>%
  filter(quarter == 4)
#> mutate: new variable 'tot_cod' (double) with 7,870 unique values and 0% NA
#> filter: removed 5,271 rows (59%), 3,595 rows remaining

spde <- make_mesh(cpue3, xy_cols = c("X", "Y"),
                  n_knots = 200, 
                  type = "kmeans", seed = 42)

cod_ind <- sdmTMB(tot_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                    temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                  data = cpue3,
                  mesh = spde,
                  family = tweedie(link = "log"),
                  spatiotemporal = "AR1",
                  spatial = "on",
                  time = "year",
                  reml = FALSE,
                  control = sdmTMBcontrol(newton_loops = 1))

fle_ind <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                    temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                  data = cpue3,
                  mesh = spde,
                  family = tweedie(link = "log"),
                  spatiotemporal = "AR1",
                  spatial = "on",
                  time = "year",
                  reml = FALSE,
                  control = sdmTMBcontrol(newton_loops = 1))

sanity(cod_ind)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
sanity(fle_ind)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large

# cod_ind2 <- sdmTMB(tot_cod ~ 0 + fyear + s(depth_sc),
#                   data = cpue3,
#                   mesh = spde,
#                   family = tweedie(link = "log"),
#                   spatiotemporal = "AR1",
#                   spatial = "on",
#                   time = "year",
#                   reml = FALSE,
#                   control = sdmTMBcontrol(newton_loops = 1))
# 
# fle_ind2 <- sdmTMB(flounder ~ 0 + fyear + s(depth_sc),
#                   data = cpue3,
#                   mesh = spde,
#                   family = tweedie(link = "log"),
#                   spatiotemporal = "AR1",
#                   spatial = "on",
#                   time = "year",
#                   reml = FALSE,
#                   control = sdmTMBcontrol(newton_loops = 1))
# 
# sanity(cod_ind2)
# sanity(fle_ind2)

Now try old data as well

# d_old <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cpue_q_1_4.csv")
# 
# # Filter to match the data I want to predict on
# d_old <- d_old %>%
#   filter(year > 1992 & year < 2020 & quarter == 4)
# 
# # Calculate standardized variables
# d_old <- d_old %>% 
#   mutate(depth_sc = (depth - mean(depth))/sd(depth),
#          X = X/1000,
#          Y = Y/1000,
#          year = as.integer(year),
#          fyear  = as.factor(year)) %>% 
#   drop_na(depth) %>% 
#   rename("density_cod" = "density") # to fit better with how flounder is named
# 
# nrow(d_old)
# nrow(cpue3)
# 
# # Plot comparison
# p1 <- ggplot() + 
#   geom_histogram(data = d_old, aes(density_cod, fill = "old"), alpha = 0.5) +
#   geom_histogram(data = cpue3, aes(tot_cod, fill = "new"), alpha = 0.5) + 
#   scale_y_continuous(trans = "log")
# 
# p2 <- ggplot() + 
#   geom_histogram(data = d_old, aes(density_fle, fill = "old"), alpha = 0.5) +
#   geom_histogram(data = cpue3, aes(flounder, fill = "new"), alpha = 0.5) + 
#   scale_y_continuous(trans = "log")
# 
# p1 + p2

# spde_old <- make_mesh(d_old, xy_cols = c("X", "Y"),
#                   n_knots = 200, 
#                   type = "kmeans", seed = 42)
# 
# cod_ind3 <- sdmTMB(density_cod ~ 0 + fyear + s(depth_sc),
#                   data = d_old,
#                   mesh = spde_old,
#                   family = tweedie(link = "log"),
#                   spatiotemporal = "AR1",
#                   spatial = "on",
#                   time = "year",
#                   reml = FALSE,
#                   control = sdmTMBcontrol(newton_loops = 1))
# 
# fle_ind3 <- sdmTMB(density_fle ~ 0 + fyear + s(depth_sc),
#                   data = d_old,
#                   mesh = spde_old,
#                   family = tweedie(link = "log"),
#                   spatiotemporal = "AR1",
#                   spatial = "on",
#                   time = "year",
#                   reml = FALSE,
#                   control = sdmTMBcontrol(newton_loops = 1))
# 
# sanity(cod_ind3)
# sanity(fle_ind3)

### Cod
# SD 24
preds_cod24_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 24), nsim = 100)
#> filter: removed 513,744 rows (87%), 79,016 rows remaining

# SD 25
preds_cod25_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 25), nsim = 100)
#> filter: removed 446,040 rows (75%), 146,720 rows remaining

# SD 26
preds_cod26_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 26), nsim = 100)
#> filter: removed 457,856 rows (77%), 134,904 rows remaining

# SD 27
preds_cod27_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 27), nsim = 100)
#> filter: removed 498,120 rows (84%), 94,640 rows remaining

# SD 28
preds_cod28_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 28), nsim = 100)
#> filter: removed 455,280 rows (77%), 137,480 rows remaining

# Now calculate the CPUE index (average)
index24_cod_sim <- get_index_sims(preds_cod24_sim, area = rep(4*4, nrow(preds_cod24_sim)))
index25_cod_sim <- get_index_sims(preds_cod25_sim, area = rep(4*4, nrow(preds_cod25_sim)))
index26_cod_sim <- get_index_sims(preds_cod26_sim, area = rep(4*4, nrow(preds_cod26_sim)))
index27_cod_sim <- get_index_sims(preds_cod27_sim, area = rep(4*4, nrow(preds_cod27_sim)))
index28_cod_sim <- get_index_sims(preds_cod28_sim, area = rep(4*4, nrow(preds_cod28_sim)))

# Add some columns
index24_cod_sim <- index24_cod_sim %>% mutate(sub_div = "24", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index25_cod_sim <- index25_cod_sim %>% mutate(sub_div = "25", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index26_cod_sim <- index26_cod_sim %>% mutate(sub_div = "26", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index27_cod_sim <- index27_cod_sim %>% mutate(sub_div = "27", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index28_cod_sim <- index28_cod_sim %>% mutate(sub_div = "28", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA


### Flounder
# SD 24
preds_fle24_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 24), nsim = 100)
#> filter: removed 513,744 rows (87%), 79,016 rows remaining

# SD 25
preds_fle25_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 25), nsim = 100)
#> filter: removed 446,040 rows (75%), 146,720 rows remaining

# SD 26
preds_fle26_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 26), nsim = 100)
#> filter: removed 457,856 rows (77%), 134,904 rows remaining

# SD 27
preds_fle27_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 27), nsim = 100)
#> filter: removed 498,120 rows (84%), 94,640 rows remaining

# SD 28
preds_fle28_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 28), nsim = 100)
#> filter: removed 455,280 rows (77%), 137,480 rows remaining

# Now calculate the CPUE index (average)
index24_fle_sim <- get_index_sims(preds_fle24_sim, area = rep(4*4, nrow(preds_fle24_sim)))
index25_fle_sim <- get_index_sims(preds_fle25_sim, area = rep(4*4, nrow(preds_fle25_sim)))
index26_fle_sim <- get_index_sims(preds_fle26_sim, area = rep(4*4, nrow(preds_fle26_sim)))
index27_fle_sim <- get_index_sims(preds_fle27_sim, area = rep(4*4, nrow(preds_fle27_sim)))
index28_fle_sim <- get_index_sims(preds_fle28_sim, area = rep(4*4, nrow(preds_fle28_sim)))

# Add some columns
index24_fle_sim <- index24_fle_sim %>% mutate(sub_div = "24", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index25_fle_sim <- index25_fle_sim %>% mutate(sub_div = "25", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index26_fle_sim <- index26_fle_sim %>% mutate(sub_div = "26", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index27_fle_sim <- index27_fle_sim %>% mutate(sub_div = "27", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index28_fle_sim <- index28_fle_sim %>% mutate(sub_div = "28", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA

div_index_sim <- bind_rows(index24_cod_sim,
                           index25_cod_sim,
                           index26_cod_sim,
                           index27_cod_sim,
                           index28_cod_sim,
                           index24_fle_sim,
                           index25_fle_sim,
                           index26_fle_sim,
                           index27_fle_sim,
                           index28_fle_sim) %>% 
  mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000) # convert to tonnes 
#> mutate: new variable 'est_t' (double) with 280 unique values and 0% NA
#>         new variable 'lwr_t' (double) with 280 unique values and 0% NA
#>         new variable 'upr_t' (double) with 280 unique values and 0% NA

Plot the index and save as csv

# # First compare with previous index:
# cod_fle_index_old <- read.csv("/Users/maxlindmark/Dropbox/Max work/R/cod_condition/output/cod_fle_index.csv") %>% 
#   mutate(sub_div = as.character(sub_div)) %>% 
#   filter(sub_div %in% c("25", "28"))
# 
# index_comp <- bind_rows(cod_fle_index_old %>% mutate(source = "cod_condition"),
#                         div_index_sim %>% mutate(source = "cod_interactions"))
# 
# ggplot(index_comp, aes(year, est_t, color = source, fill = source)) +
#   facet_wrap(sub_div ~ species, scales = "free") +
#   geom_line() +
#   geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t), alpha = 0.2, color = NA)
# 
# 
# ## If it isn't the same with the data, plot pred grids against each other

ggplot(div_index_sim, aes(year, est_t, color = sub_div)) +
  geom_line() +
  facet_wrap(~species, scales = "free") +
  geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2, color = NA) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
  guides(color = "none") + 
  labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
  NULL

div_index_sim %>% 
  filter(sub_div %in% c(25, 28)) %>% 
  ggplot(aes(year, est_t/1000, color = species, fill = species)) +
  geom_line() +
  facet_wrap(~sub_div, scales = "free", ncol = 1) +
  geom_ribbon(aes(year, ymin = lwr_t/1000, ymax = upr_t/1000), alpha = 0.2, color = NA) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "Species") +
  guides(color = "none") + 
  labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
  NULL
#> filter: removed 168 rows (60%), 112 rows remaining

ggplot(div_index_sim, aes(year, est_t, color = sub_div)) +
  geom_line() +
  facet_wrap(~species, scales = "free") +
  #geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
  guides(color = "none") + 
  labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
  NULL

ggplot(div_index_sim, aes(year, est_t, color = species)) +
  geom_line() +
  facet_wrap(~sub_div, scales = "free") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
  labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
  NULL

write.csv(div_index_sim, "output/cod_fle_index.csv")